
## PEASANT RESISTANCE AND ECONOMIC AFFLUENCE: LESSONS FROM PARAGUAY
## Appendix A: Additional Figures and Tables
## January 2024

## Code Authors:
## Germán Feierherd (gfeierherd@udesa.edu.ar) 
## Jorge Mangonnet (jorge.mangonnet@vanderbilt.edu)


#### Preamble ------------------------------------------------------------------

## Workspace setup
rm(list = ls())		         # clear all objects in memory
options(max.print = 5000)  # expand display
options(scipen = 999)      # disable sci notation
dev.off()                  # reload graphic device
cat("\014")                # clear console

## Load/install packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  readxl,
  interflex,
  marginaleffects,
  haven,
  readr,
  reshape2,
  magrittr,
  dplyr,
  plm,
  lfe,
  pcse,
  lmtest,
  stargazer,
  xtable,
  texreg,
  estimatr, 
  fixest,
  margins, 
  sjPlot, 
  multiwayvcov,
  #  rgdal,
  ggplot2,
  wesanderson, 
  fixest,
  did,
  sf)

## Set working directory
gf <- "~/Dropbox/Working_Papers/investigacion protesta Paraguay/"; main <- gf  # german's path
jm <- "~/Dropbox/Research/investigacion protesta Paraguay/"; main <- jm        # jorge's path
setwd(main)

## Functions
source(paste0(main, "4_code/functions_paraguay.R"))

## Set CRS for shapefiles 
wgs84 <- CRS("+proj=longlat +datum=WGS84 +no_defs")

## Load full dataset
load(paste0(main, "4_code/DataSetPY.Rda"))


#### Appendix A: Additional Figures and Tables ---------------------------------

#### Figure A.2: Prices of Capital-Intensive Crops, 1991-2017 ----

# (a) Soybeans
price.soybean <- read.csv(paste0(main, "/3_sources/fao/FAOSTAT_data_Soybean_Prices.csv")) %>% 
  dplyr::select(Year, Value) %>%
  dplyr::rename(year = Year, soybean_price = Value)

price.soybean <- ggplot(price.soybean, aes(year, soybean_price)) +
  geom_line() +
  geom_line(data = subset(price.soybean, year > 1999 & year < 2015), color = 'red') +
  theme_bw() + 
  coord_cartesian(ylim = c(5,460)) + 
  theme(panel.border = element_blank(),
        axis.text.y = element_text(size = 14, hjust = 1),
        axis.text.x = element_text(size = 14, angle = 45),
        axis.title = element_text(size = 16, hjust =.5)) + 
  geom_vline(xintercept = 2000, linetype = "dashed", color = "red") +
  geom_vline(xintercept = 2014, linetype = "dashed", color = "red") +
  scale_x_continuous(breaks = (seq(1991, 2017, 1))) + 
  xlab("Year") + ylab("Soybean Prices (USD per metric ton)")

pdf(file = paste0(main, "1_writing/2023/figures/price_soybean.pdf"))
price.soybean
dev.off()

# (b) Maize
price.maize <- read.csv(paste0(main, "/3_sources/fao/FAOSTAT_data_Maize_Prices.csv")) %>% 
  dplyr::select(Year, Value) %>%
  dplyr::rename(year = Year, maize_price = Value)

price.maize <- ggplot(price.maize, aes(year, maize_price)) + 
  geom_line() +
  geom_line(data = subset(price.maize, year > 1999 & year < 2015), color = 'red') +
  theme_bw() + 
  coord_cartesian(ylim = c(5,460)) + 
  theme(panel.border = element_blank(),
        axis.text.y = element_text(size = 14, hjust = 1),
        axis.text.x = element_text(size = 14, angle = 45),
        axis.title = element_text(size = 16, hjust =.5)) + 
  geom_vline(xintercept = 2000, linetype = "dashed", color = "red") +
  geom_vline(xintercept = 2014, linetype = "dashed", color = "red") +
  scale_x_continuous(breaks = (seq(1991, 2017, 1))) + 
  xlab("Year") + ylab("Maize Prices (USD per metric ton)")

pdf(file = paste0(main, "1_writing/2023/figures/price_maize.pdf"))
price.maize
dev.off()

# (c) Sugar
price.sugar <- read.csv(paste0(main, "/3_sources/fao/FAOSTAT_data_Sugar_Prices.csv")) %>% 
  dplyr::select(Year, Value) %>%
  dplyr::rename(year = Year, sugar_price = Value)

price.sugar <- ggplot(price.sugar, aes(year, sugar_price)) + 
  geom_line() +
  geom_line(data = subset(price.sugar, year > 1999 & year < 2015), color = 'red') +
  theme_bw() + 
  coord_cartesian(ylim = c(5,460)) + 
  theme(panel.border = element_blank(),
        axis.text.y = element_text(size = 14, hjust = 1),
        axis.text.x = element_text(size = 14, angle = 45),
        axis.title = element_text(size = 16, hjust =.5)) + 
  geom_vline(xintercept = 2000, linetype = "dashed", color = "red") +
  geom_vline(xintercept = 2014, linetype = "dashed", color = "red") +
  scale_x_continuous(breaks = (seq(1991, 2017, 1))) + 
  xlab("Year") + ylab("Sugar Prices (USD per metric ton)")

pdf(file = paste0(main, "1_writing/2023/figures/price_sugar.pdf"))
price.sugar
dev.off()

remove(price.comm, price.soybean, price.maize, price.rice, price.sugar, price.cotton)


#### Figure A.3: Local Regression of Peasant Resistance on Land Suitability ----

# 2000
lwr.2000 <- ggplot(dta[which(dta$year == 2000), ], aes(ln_suitability, ln_conflict)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") + 
  coord_cartesian(ylim = c(-.2,1)) +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 12, hjust = 1),
        axis.title = element_text(size = 15, hjust = .5)) +
  ggtitle("2000") + 
  xlab(" ") + ylab(" ")

# 2001
lwr.2001 <- ggplot(dta[which(dta$year == 2001), ], aes(ln_suitability, ln_conflict)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") + 
  coord_cartesian(ylim = c(-.2,1)) +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 12, hjust = 1),
        axis.title = element_text(size = 15, hjust = .5)) +
  ggtitle("2001") + 
  xlab(" ") + ylab(" ")

# 2002
lwr.2002 <- ggplot(dta[which(dta$year == 2002), ], aes(ln_suitability, ln_conflict)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") + 
  coord_cartesian(ylim =c(-.2,1)) +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 12, hjust = 1),
        axis.title = element_text(size = 15, hjust = .5)) +
  ggtitle("2002") + 
  xlab(" ") + ylab(" ")

# 2003
lwr.2003 <- ggplot(dta[which(dta$year == 2003), ], aes(ln_suitability, ln_conflict)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") + 
  coord_cartesian(ylim = c(-.2,1)) +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 12, hjust = 1),
        axis.title = element_text(size = 15, hjust = .5)) +
  ggtitle("2003") + 
  xlab(" ") + ylab(" ")

# 2004
lwr.2004 <- ggplot(dta[which(dta$year == 2004), ], aes(ln_suitability, ln_conflict)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") + 
  coord_cartesian(ylim=c(-.2,1)) +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 12, hjust = 1),
        axis.title = element_text(size = 15, hjust = .5)) +
  ggtitle("2004") + 
  xlab(" ") + ylab(" ")

# 2005
lwr.2005 <- ggplot(dta[which(dta$year == 2005), ], aes(ln_suitability, ln_conflict)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") + 
  coord_cartesian(ylim = c(-.2,1)) +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 12, hjust = 1),
        axis.title = element_text(size = 15, hjust = .5)) +
  ggtitle("2005") + 
  xlab(" ") + ylab(" ")

# 2006 -- add y-axis label
lwr.2006 <- ggplot(dta[which(dta$year == 2006), ], aes(ln_suitability, ln_conflict)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") + 
  coord_cartesian(ylim = c(-.2,1)) +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 11, hjust = 1),
        axis.title = element_text(size = 11, hjust = .5)) +
  ggtitle("2006") + 
  xlab(" ") + ylab("Events of\nPeasant Resistance")

# 2007
lwr.2007 <- ggplot(dta[which(dta$year == 2007), ], aes(ln_suitability, ln_conflict)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha = .2), sides = "b") + 
  coord_cartesian(ylim = c(-.2,1)) +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 12, hjust = 1),
        axis.title = element_text(size = 15, hjust = .5)) +
  ggtitle("2007") + 
  xlab(" ") + ylab(" ")

# 2008
lwr.2008 <- ggplot(dta[which(dta$year == 2008), ], aes(ln_suitability, ln_conflict)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha = .2), sides = "b") + 
  coord_cartesian(ylim = c(-.2,1)) +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 12, hjust = 1),
        axis.title = element_text(size = 15, hjust = .5)) +
  ggtitle("2008") + 
  xlab(" ") + ylab(" ")

# 2009
lwr.2009 <- ggplot(dta[which(dta$year == 2009), ], aes(ln_suitability, ln_conflict)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha = .2), sides = "b") + 
  coord_cartesian(ylim = c(-.2,1)) +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 12, hjust = 1),
        axis.title = element_text(size = 15, hjust = .5)) +
  ggtitle("2009") + 
  xlab(" ") + ylab(" ")

# 2010
lwr.2010 <- ggplot(dta[which(dta$year == 2010), ], aes(ln_suitability, ln_conflict)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha = .2), sides = "b") + 
  coord_cartesian(ylim = c(-.2,1)) +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 12, hjust = 1),
        axis.title = element_text(size = 15, hjust = .5)) +
  ggtitle("2010") + 
  xlab(" ") + ylab(" ")

# 2011
lwr.2011 <- ggplot(dta[which(dta$year == 2011), ], aes(ln_suitability, ln_conflict)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha = .2), sides = "b") + 
  coord_cartesian(ylim = c(-.2,1)) +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 12, hjust = 1),
        axis.title = element_text(size = 15, hjust = .5)) +
  ggtitle("2011") + 
  xlab(" ") + ylab(" ")

# 2012
lwr.2012 <- ggplot(dta[which(dta$year == 2012), ], aes(ln_suitability, ln_conflict)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") + 
  coord_cartesian(ylim = c(-.2,1)) +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 12, hjust = 1),
        axis.title = element_text(size = 15, hjust = .5)) +
  ggtitle("2012") + 
  xlab(" ") + ylab(" ")

# 2013 
lwr.2013 <- ggplot(dta[which(dta$year == 2013), ], aes(ln_suitability, ln_conflict)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") + 
  coord_cartesian(ylim = c(-.2,1)) +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 11, hjust = 1),
        axis.title = element_text(size = 11, hjust = .5)) +
  ggtitle("2013") + 
  xlab("Land Suitability") + ylab(" ")

# 2014 
lwr.2014 <- ggplot(dta[which(dta$year == 2014), ], aes(ln_suitability, ln_conflict)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") + 
  coord_cartesian(ylim = c(-.2,1)) +
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 11, hjust = 1),
        axis.title = element_text(size = 11, hjust = .5)) +
  ggtitle("2014") + 
  xlab(" ") + ylab(" ")

# Arrange plots
lwr.plots <- ggarrange(
  lwr.2000, lwr.2001, lwr.2002, 
  lwr.2003, lwr.2004, lwr.2005, 
  lwr.2006, lwr.2007, lwr.2008, 
  lwr.2009, lwr.2010, lwr.2011, 
  lwr.2012, lwr.2013, lwr.2014, 
  ncol = 3, nrow = 5)

pdf(file = paste0(main, "1_writing/2023/figures/loess_plots.pdf"))
lwr.plots
dev.off()

remove(lwr.2000, lwr.2001, lwr.2002, 
       lwr.2003, lwr.2004, lwr.2005, 
       lwr.2006, lwr.2007, lwr.2008, 
       lwr.2009, lwr.2010, lwr.2011, 
       lwr.2012, lwr.2013, lwr.2014,
       lwr.plots)


#### Figure A.5: Subsistence Settlements (núcleos) by Municipality, 1990-1992 ----
#### Figure A.6: Peasant Committees by Municipality, 1992-1993 ----

# Group the data
sa <- dta %>%
  dplyr::select(c(codigo, subsistence_farms_median)) %>% 
  group_by(codigo) %>%
  summarise(subsistence_farms_median = median(subsistence_farms_median)) %>%
  mutate(subsistence_farms_median = as.factor(subsistence_farms_median))

po <- dta %>%
  dplyr::select(c(codigo, orgs_locales_median)) %>% 
  group_by(codigo) %>%
  summarise(orgs_locales_median = median(orgs_locales_median)) %>%
  mutate(orgs_locales_median = replace_na(orgs_locales_median, 0),
         orgs_locales_median = as.factor(orgs_locales_median))

# Load municipal polygons and merge
py <- readOGR(paste0(main, "7_gis/paraguay_2002_municipios"), "paraguay_2002_municipios") %>%
  spTransform(wgs84) %>%
  subset(select = codigo)
py@data$id <- seq(1, 224, 1)                    # ids
py@data$codigo <- as.character(py@data$codigo)  # codigos

# Merge and get coordinates
py <- merge(x = py, y = sa, by = "codigo") %>%
  broom::tidy(region = "codigo")
pds <- fortify(py, region = "codigo")


## (a) Susbsistence settlements (Figure A.5) --

map.subsistence <- ggplot() + theme_void() + coord_equal() + 
  geom_map(data = py, aes(x = long, y = lat, map_id = id), map = pds, colour = "black", fill = NA, size = 1) +
  geom_map(data = sa, map = py, aes(fill = subsistence_farms_median, map_id = codigo, color = subsistence_farms_median)) + 
  scale_fill_manual(name = "Subsistence\nAgriculture", labels = c("Low", "High"), values = c("lightgrey", "green")) + 
  scale_color_manual(guide = FALSE, values = c("black", "black")) +
  theme(legend.position = "right") + 
  theme(legend.text = element_text(size = 14)) +
  theme(legend.title = element_text(size = 14)) + 
  #ggtitle("Subsistence Agriculture") +
  #theme(plot.title = element_text(size = 18)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))

png(file = paste0(main, "1_writing/2023/figures/map_subsist.png"), width = 8, height = 8, res = 300, units='in')
map.subsistence
dev.off()

closeAllConnections()


## (b) Peasant committees (Figure A.6) --

map.organization <- ggplot() + theme_void() + coord_equal() + 
  geom_map(data = py, aes(x = long, y = lat, map_id = id), map = pds, colour = "black", fill = NA, size = 1) +
  geom_map(data = po, map = py, aes(fill = orgs_locales_median, map_id = codigo, color = orgs_locales_median)) + 
  scale_fill_manual(name = "Organizational\nCapacities", labels = c("Low", "High"), values = c("lightgrey", "red")) + 
  scale_color_manual(guide = FALSE, values = c("black", "black")) +
  theme(legend.position = "right") + 
  theme(legend.text = element_text(size = 14)) +
  theme(legend.title = element_text(size = 14)) + 
  #ggtitle("Subsistence Agriculture") +
  #theme(plot.title = element_text(size = 18)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))

png(file = paste0(main, "1_writing/2023/figures/map_organiz.png"), width = 8, height = 8, res = 300, units='in')
map.organization
dev.off()

closeAllConnections()

rm(map.organization,map.subsistence)


#### Table A.1: Descriptive Statistics ----

# Select variables
desc.vars <- subset(dta, select = c(
  ln_conflict,
  b_conflict, ln_positive,
  ln_conflict_t1, ln_conflict_t2, ln_conflict_t3,
  n_conflict,
  ln_nochaco,
  ln_conflict_inc,
  ln_peasvsindig,
  ln_landvsstate,
  ln_all,
  ln_price, ln_suitability,
  ln_price_pbo, ln_suitability_pbo,
  inflation,
  ln_soy_price, ln_mze_price, ln_sug_price,
  ln_soy_suit, ln_mze_suit, ln_sug_suit,
  subsistence_farms_median, orgs_locales_median,
  subsistence_farms_1Q, orgs_locales_1Q,
  subsistence_farms_3Q, orgs_locales_3Q,
  cotton_farms20_median,
  old_colonies_median,
  landowner_orgs_median,
  total_planted_diff,
  soybean_planted_diff, maize_planted_diff, sugar_planted_diff, 
  soybean_planted_diff_14,
  distance_hubs, distance_asuncion, untitled_share, ln_pop))

# Define labels  
desc.labels <- c("Resistance",
                 "Resistance (binary)",
                 "Resistance (restricted sample)",
                 "Resistance (t-1)", 
                 "Resistance (t-2)", 
                 "Resistance (t-3)",
                 "Resistance (discrete)",
                 "Resistance (w/o Chaco)",
                 "Resistance (income events)",
                 "Resistance (peasant vs. indigenous)",
                 "Landowner Protest",
                 "All Rural Unrest",
                 "Prices", 
                 "Suitability",
                 "Prices (labor-intensive)", 
                 "Suitability (labor-intensive)",
                 "Consumer Price Index",
                 "Soybean Prices", 
                 "Maize Prices", 
                 "Sugar Prices",
                 "Soybean Suitability", 
                 "Maize Suitability", 
                 "Sugar Suitability", 
                 "Subsistence Settlements", "Peasant Committees",
                 "Subsistence Settlements (1Q)", "Subsistence Settlements (3Q)",
                 "Peasant Committees (1Q)", "Peasant Committees (3Q)",
                 "Peasant Cotton Farms",
                 "IBR Colonies",
                 "Landowner Associations",
                 #"Evictions (procedures)", "Evictions (families)", "Evictions (hectares)",
                 "Total Planted (2008-1991)",
                 "Soybean Planted (2008-1991)", "Maize Planted (2008-1991)", "Sugar Planted (2008-1991)", 
                 "Soybean Planted (2014-2008)", 
                 "Distance from Export Hubs", "Distance from Asunci\\'{o}n", "Untitled Farmland", "Population"
)

# Design table
names(desc.vars) <- desc.labels
desc.stats.table <- psych::describe(desc.vars, fast = TRUE, quant = c(.25, .5, .75))
desc.stats.table[[1]] <- NULL
desc.stats.table[[6]] <- NULL
desc.stats.table[[6]] <- NULL
desc.stats.table[[6]] <- NULL
desc.stats.table[[6]] <- NULL
desc.stats.table[[6]] <- NULL
colnames(desc.stats.table)[1] <- "N"
colnames(desc.stats.table)[2] <- "Mean"
colnames(desc.stats.table)[3] <- "Std.Dev"
colnames(desc.stats.table)[4] <- "Min"
colnames(desc.stats.table)[5] <- "Max"
desc.stats.table[[1]] <- round(desc.stats.table[[1]], 0)
desc.stats.table$N <- as.character(desc.stats.table$N)

sink(paste0(main,"1_writing/2023/tables/table_descstats.tex"))
print(xtable(desc.stats.table, type = "latex"), size = "small", floating = FALSE)
sink()

closeAllConnections()


#### Figure A.6: Event Study ----

# Create time indicator
dta <- dta %>%
  group_by(codigo) %>%
  mutate(time = ifelse(year %in% c(2000), -3, 
                       ifelse(year %in% c(2001, 2005), -2,
                              ifelse(year %in% c(2002, 2006, 2010), -1, 
                                     ifelse(year %in% c(2003, 2007, 2011), 0,
                                            ifelse(year %in% c(2004, 2008, 2012), 1, 
                                                   ifelse(year %in% c(2013), 2, 3))))))) # Small typo here I think, changed 2003 for 2004

# Event study model
es <- feols(ln_conflict ~ i(time, ln_suitability, -1) | codigo + time, cluster = "codigo", data = dta)

# Dataset
dta_es <- data.frame(
  b = es[["coefficients"]], 
  se = es[["se"]], 
  ci.u = es[["coefficients"]] + (1.96*es[["se"]]), 
  ci.l = es[["coefficients"]] - (1.96*es[["se"]])) %>%
  dplyr::add_row(b = 0, se = 0, ci.u = 0, ci.l = 0, .before = 3) %>%
  mutate(time = c(-3, -2, -1, 0, 1, 2, 3))

# Plot
pdf(file = paste0(main, "1_writing/2023/figures/event_study_2014.pdf"), width = 5, height = 5)

plot(dta_es$time, dta_es$b, 
     type = "p", pch = 16, lwd = 1, 
     cex = .75, cex.lab = 1, axes = TRUE, ylim = c(-2.75,.75), 
     xlab = "Years Before/After Price Hike", 
     ylab = "Effect of Land Suitability",
     frame.plot = FALSE)

arrows(dta_es$time, dta_es$ci.l, dta_es$time, dta_es$ci.u, length = 0.025, angle = 90, code = 3)

abline(v = -1, lty = 2, col="red")
abline(h = 0, lty = 3)

dev.off()

